qsmooth test

library(qsmooth)

Setting options('download.file.method.GEOquery'='auto')
Setting options('GEOquery.inmemory.gpl'=FALSE)

library(Rtsne)
library(ggplot2)
library(formattable)
library(tidyverse)

Loading tidyverse: tibble
Loading tidyverse: tidyr
Loading tidyverse: readr
Loading tidyverse: purrr
Loading tidyverse: dplyr
Conflicts with tidy packages -------------------------------------------------------------------------------------------------------------------------------------------------------------------
filter(): dplyr, stats
lag():    dplyr, stats

source('~/git/scripts/theme_Publication.R')
# https://github.com/stephaniehicks/qsmooth
load('~/git/unified_gene_expression/data/lengthScaledTPM_eye_gtex.Rdata')
source('~/git/unified_gene_expression/scripts/parse_sample_attribute.R')

data.table 1.9.6  For help type ?data.table or https://github.com/Rdatatable/data.table/wiki
The fastest way to learn (by data.table authors): https://www.datacamp.com/courses/data-analysis-the-data-table-way
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
data.table + dplyr code now lives in dtplyr.
Please library(dtplyr)!
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Attaching package: ‘data.table’

The following objects are masked from ‘package:dplyr’:

    between, last

The following object is masked from ‘package:purrr’:

    transpose

Testing to see if running qsmooth (https://github.com/stephaniehicks/qsmooth, weighted quantile normalization) improves the clustering performance.

lengthScaledTPM <- lengthScaledTPM[,!(is.na(lengthScaledTPM[1,]))]
samples <- data.frame(colnames(lengthScaledTPM))
colnames(samples) <- 'sample_accession'
tissue_frame <- left_join(samples, core_info) %>% select(sample_accession, Tissue, Sub_Tissue) %>% distinct()

Joining, by = "sample_accession"
joining character vector and factor, coercing into character vector

# make sure they are lined up
tissues <- left_join(samples, tissue_frame)

Joining, by = "sample_accession"
joining character vector and factor, coercing into character vector

# just keep the keeper set of tissues (no prostrate, etc)
sub <- tissues %>% filter(Tissue %in% keepers)
lengthScaledTPM_sub<-lengthScaledTPM[,sub$sample_accession]
qs <- qsmooth(object = lengthScaledTPM_sub,groupFactor = as.factor(sub$Tissue))
lengthScaledTPM_qsmooth <- qsmoothData(qs)
set.seed(23235)
tsne_out <- Rtsne(as.matrix(log2(t(lengthScaledTPM_qsmooth)+1)),perplexity = 40, check_duplicates = FALSE, theta=0.0 )
tsne_plot <- data.frame(tsne_out$Y)
tsne_plot$sample_accession <- colnames(lengthScaledTPM_qsmooth)
tsne_plot %>% left_join(.,core_info)  %>%
  ggplot(.,aes(x=X1,y=X2,colour=Tissue,shape=Tissue)) + 
  geom_point(size=4) + scale_shape_manual(values=c(0:20,35:50)) +
  ggtitle(paste0("t-sne. Perplexity = ", 40)) +
  theme_Publication()

Joining, by = "sample_accession"

No qsmooth

set.seed(23235)
tsne_out <- Rtsne(as.matrix(log2(t(lengthScaledTPM_sub)+1)),perplexity = 40, check_duplicates = FALSE, theta=0.0 )
tsne_plot <- data.frame(tsne_out$Y)
tsne_plot$sample_accession <- colnames(lengthScaledTPM_sub)
tsne_plot %>% left_join(.,core_info)  %>%
  ggplot(.,aes(x=X1,y=X2,colour=Tissue,shape=Tissue)) + 
  geom_point(size=4) + scale_shape_manual(values=c(0:20,35:50)) +
  ggtitle(paste0("t-sne. Perplexity = ", 40)) +
  theme_Publication()

Joining, by = "sample_accession"

Expression profiles without qsmooth

gather_lST<-gather(lengthScaledTPM_sub,sample_accession) %>% left_join(.,core_tight)

Joining, by = "sample_accession"

ggplot(gather_lST,aes(x=log2(value+1),group=sample_accession))+geom_density()+facet_wrap(~Tissue)+coord_cartesian(ylim=c(0,0.5)) + theme_Publication()

Expression profiles with qsmooth. Odd spikes in the Cornea and RPE. Some weird patterns in number of lower-expressed genes (see Blood Vessel, Brain, Heart)

core_tight2 <- core_tight
core_tight2$sample_accession <- gsub(pattern='E-MTAB-',replacement = 'E.MTAB.',core_tight$sample_accession)
gather_lST<-gather(data.frame(lengthScaledTPM_qsmooth),sample_accession) %>% left_join(.,core_tight2)

Joining, by = "sample_accession"

ggplot(gather_lST,aes(x=log2(value+1),group=sample_accession))+geom_density()+facet_wrap(~Tissue)+coord_cartesian(ylim=c(0,0.5))+theme_Publication()

Let’s look at just the Cornea. Funky.

core_tight2 <- core_tight
core_tight2$sample_accession <- gsub(pattern='E-MTAB-',replacement = 'E.MTAB.',core_tight$sample_accession)
gather_lST <- gather(data.frame(lengthScaledTPM_qsmooth),sample_accession) %>% left_join(.,core_tight2) %>% filter(Tissue=='Cornea')

Joining, by = "sample_accession"

ggplot(gather_lST,aes(x=log2(value+1),group=sample_accession))+geom_density()+facet_wrap(~Tissue)+coord_cartesian(ylim=c(0,0.5)) + theme_Publication()

Let’s find the weird ones. SRS390607 to SRS390610. Weird median. Plotting the uncorrected values…Huh, seems like there are lots of zeros….

gather(data.frame(lengthScaledTPM_qsmooth),sample_accession) %>% left_join(.,core_tight2) %>% filter(Tissue=='Cornea') %>% group_by(sample_accession) %>% summarise(median(value)) %>% arrange(-`median(value)`) %>% head() %>% formattable()

Joining, by = "sample_accession"

core_eye_info %>% filter(sample_accession=='SRS390607') %>% formattable()

gather_lST<-gather(lengthScaledTPM_sub,sample_accession) %>% left_join(.,core_tight) %>% filter(sample_accession %in% c('SRS390607','SRS390608','SRS390609','SRS390610'))

Joining, by = "sample_accession"

ggplot(gather_lST,aes(x=log2(value+1)))+geom_density()+facet_wrap(~Tissue)+coord_cartesian(ylim=c(0,0.5)) + theme_Publication()

Let’s look at the number level for one of those odd samples. Uh, yikes. Let’s check the median for ALL of the samples to see if any others have this problem. I’m guessing more do.

summary(lengthScaledTPM_sub[,'SRS390607'])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       0     320       0 6503000 

Yeah, a few. A plateau at a median of 50.

apply(lengthScaledTPM,2,function(x) median(x)) %>% sort() %>% head(60) %>% data.frame() %>% formattable()

apply(lengthScaledTPM,2,function(x) median(x)) %>% density() %>% plot()
axis(side = 1, at=seq(0,500,10))

Let’s redo the qsmooth’ed density plots above, but with the low (<10) colored. Bam, that’s it.

median_TPM <- apply(lengthScaledTPM,2,function(x) median(x)) %>% sort() %>% formattable()
lows <- names(median_TPM[median_TPM<10])
core_tight2 <- core_tight
core_tight2$sample_accession <- gsub(pattern='E-MTAB-',replacement = 'E.MTAB.',core_tight$sample_accession)
gather_lST<-gather(data.frame(lengthScaledTPM_qsmooth),sample_accession) %>% left_join(.,core_tight2) %>% mutate(LowMedian=ifelse(sample_accession %in% lows, 'LowMedian','OK'))

Joining, by = "sample_accession"

ggplot(gather_lST,aes(x=log2(value+1),group=sample_accession,colour=LowMedian))+geom_density(alpha=0.5)+facet_wrap(~Tissue)+coord_cartesian(ylim=c(0,0.5))+theme_Publication()

Let’s up the median threshold to, say, 50. This is probably a reasonable place to set it.

median_TPM <- apply(lengthScaledTPM,2,function(x) median(x)) %>% sort() %>% formattable()
lows <- names(median_TPM[median_TPM<60])
core_tight2 <- core_tight
core_tight2$sample_accession <- gsub(pattern='E-MTAB-',replacement = 'E.MTAB.',core_tight$sample_accession)
gather_lST<-gather(data.frame(lengthScaledTPM_qsmooth),sample_accession) %>% left_join(.,core_tight2) %>% mutate(LowMedian=ifelse(sample_accession %in% lows, 'LowMedian','OK'))

Joining, by = "sample_accession"

ggplot(gather_lST,aes(x=log2(value+1),group=sample_accession,colour=LowMedian))+geom_density(alpha=0.5)+facet_wrap(~Tissue)+coord_cartesian(ylim=c(0,0.5))+theme_Publication()

Ugh, this tosses a LOT of my fetal RPE (44 originally). Eh, life.

gather_lST %>% filter(sample_accession %in% lows) %>% select(sample_accession, Tissue, Sub_Tissue) %>% distinct() %>% group_by(Tissue, Sub_Tissue) %>% summarise(Count=n()) %>% formattable()

Redo the tsne (qsmooth’ed) and tossing the above samples. Looks better. A couple of Retina ($) samples that have wandered away (X2=10,X1=-10).

set.seed(23235)
lengthScaledTPM_qsmooth_highExp <- as.matrix(lengthScaledTPM_qsmooth[,!colnames(lengthScaledTPM_qsmooth) %in% lows])
tsne_out <- Rtsne(as.matrix(log2(t(lengthScaledTPM_qsmooth_highExp)+1)),perplexity = 50, check_duplicates = FALSE, theta=0.0 )
tsne_plot <- data.frame(tsne_out$Y)
tsne_plot$sample_accession <- colnames(lengthScaledTPM_qsmooth_highExp)
tsne_plot %>% left_join(.,core_info)  %>%
  ggplot(.,aes(x=X1,y=X2,colour=Tissue,shape=Tissue)) + 
  geom_point(size=4) + scale_shape_manual(values=c(0:20,35:50)) +
  ggtitle(paste0("t-sne. Perplexity = ", 50)) +
  theme_Publication()

Joining, by = "sample_accession"

Let’s do a quick check on genes with overall low experession. Less than 1 count on average across all samples, for a gene.

density(rowSums(log2(lengthScaledTPM+1))) %>% plot()

table((rowSums(lengthScaledTPM)/ncol(lengthScaledTPM))<1)

FALSE  TRUE 
19277  1053 

lowly_expressed_genes <- row.names(lengthScaledTPM[(rowSums(lengthScaledTPM)/ncol(lengthScaledTPM)<1),])

OK, let’s also toss these and see how that changes the t-sne clustering.

set.seed(23235)
lengthScaledTPM_qsmooth_highExp <- as.matrix(lengthScaledTPM_qsmooth[,!colnames(lengthScaledTPM_qsmooth) %in% lows])
lengthScaledTPM_qsmooth_highExp_remove_lowGenes <- lengthScaledTPM_qsmooth_highExp[!(row.names(lengthScaledTPM_qsmooth_highExp) %in% lowly_expressed_genes),]
tsne_out <- Rtsne(as.matrix(log2(t(lengthScaledTPM_qsmooth_highExp_remove_lowGenes)+1)),perplexity = 50, check_duplicates = FALSE, theta=0.0 )
tsne_plot <- data.frame(tsne_out$Y)
tsne_plot$sample_accession <- colnames(lengthScaledTPM_qsmooth_highExp_remove_lowGenes)
tsne_plot %>% left_join(.,core_info)  %>%
  ggplot(.,aes(x=X1,y=X2,colour=Tissue,shape=Tissue)) + 
  geom_point(size=4) + scale_shape_manual(values=c(0:20,35:50)) +
  ggtitle(paste0("t-sne. Perplexity = ", 50)) +
  theme_Publication()

Joining, by = "sample_accession"

OK, this is getting hard/impossible to see whether I’m doing any better. Need to get more analytical. Metric: cluster purity. K-means cluster (k number is number of sub-tissues or tissues?), check purity of each cluster.

load('~/git/unified_gene_expression/data/lengthScaledTPM_eye_gtex.Rdata')
lengthScaledTPM <- lengthScaledTPM[,!(is.na(lengthScaledTPM[1,]))]
sub <- tissues %>% filter(Tissue %in% keepers)
lengthScaledTPM_sub<-lengthScaledTPM[,sub$sample_accession]
set.seed(23235)
tsne_out <- Rtsne(as.matrix(log2(t(lengthScaledTPM_sub)+1)),perplexity = 40, check_duplicates = FALSE, theta=0.0 )
tsne_plot <- data.frame(tsne_out$Y)
tsne_plot$sample_accession <- colnames(lengthScaledTPM_sub)
tsne_go <- tsne_plot %>% left_join(.,core_info)

Joining, by = "sample_accession"

# Determine number of clusters
wss <- (nrow(tsne_go[,c('X1','X2')])-1)*sum(apply(tsne_go[,c('X1','X2')],2,var))
for (i in 2:60) wss[i] <- sum(kmeans(tsne_go[,c('X1','X2')], 
    centers=i)$withinss)
plot(1:60, wss, type="b", xlab="Number of Clusters",
  ylab="Within groups sum of squares")

# 24 is the total number of Tissues
# 48 is the total number of Sub Tissues
set.seed(23235)
fit <- kmeans(tsne_go[,c('X1','X2')], 24)
cbind(tsne_go, fit$cluster) %>% group_by(fit$cluster,Tissue) %>% summarise(Count = n()) %>% mutate(freq = Count /sum(Count)) %>% formattable()

Now go deeper and count number of Tissue types and purity (highest freq) for each cluster

cbind(tsne_go, fit$cluster) %>% group_by(fit$cluster,Tissue) %>% summarise(Count = n()) %>% mutate(freq = Count /sum(Count)) %>% summarise(Tissue_Count = n(), highest_freq=max(freq)) %>% formattable() 

OK, now take the mean of Tissue_Count and highest_freq, for two metrics for each tnse.

cbind(tsne_go, fit$cluster) %>% group_by(fit$cluster,Tissue) %>% summarise(Count = n()) %>% mutate(freq = Count /sum(Count)) %>% summarise(Tissue_Count = n(), highest_freq=max(freq)) %>% mutate(mean(Tissue_Count), mean(highest_freq)) %>% select(`mean(Tissue_Count)`, `mean(highest_freq)`) %>% distinct()

Cool. Let’s go the end (qsmooth, removal of low-median samples, removal of low expression genes) and see if these metrics get better. YES THEY DO. IT WORKS. OMG. (Sorry, usually stuff like this blows up in my face once I test it rigorously).

set.seed(23235)
tsne_out <- Rtsne(as.matrix(log2(t(lengthScaledTPM_qsmooth_highExp_remove_lowGenes)+1)),perplexity = 40, check_duplicates = FALSE, theta=0.0 )
tsne_plot <- data.frame(tsne_out$Y)
tsne_plot$sample_accession <- colnames(lengthScaledTPM_qsmooth_highExp_remove_lowGenes)
tsne_go <- tsne_plot %>% left_join(.,core_info)

Joining, by = "sample_accession"

set.seed(23235)
fit <- kmeans(tsne_go[,c('X1','X2')], 24)
cbind(tsne_go, fit$cluster) %>% group_by(fit$cluster,Tissue) %>% summarise(Count = n()) %>% mutate(freq = Count /sum(Count)) %>% summarise(Tissue_Count = n(), highest_freq=max(freq)) %>% mutate(mean(Tissue_Count), mean(highest_freq)) %>% select(`mean(Tissue_Count)`, `mean(highest_freq)`) %>% distinct()

Let’s see whether having qsmooth correct for Sub Tissues improves performance.

set.seed(23235)
tsne_out <- Rtsne(as.matrix(log2(t(lengthScaledTPM_qsmooth_subTissue)+1)),perplexity = 40, check_duplicates = FALSE, theta=0.0 )
tsne_plot <- data.frame(tsne_out$Y)
tsne_plot$sample_accession <- colnames(lengthScaledTPM_qsmooth_subTissue)
tsne_go <- tsne_plot %>% left_join(.,core_info)

Joining, by = "sample_accession"

set.seed(23235)
fit <- kmeans(tsne_go[,c('X1','X2')], 24)
cbind(tsne_go, fit$cluster) %>% group_by(fit$cluster,Tissue) %>% summarise(Count = n()) %>% mutate(freq = Count /sum(Count)) %>% summarise(Tissue_Count = n(), highest_freq=max(freq)) %>% mutate(mean(Tissue_Count), mean(highest_freq)) %>% select(`mean(Tissue_Count)`, `mean(highest_freq)`) %>% distinct()

Now with the corrections….kind of middling. May not be enough tissues per group.

set.seed(23235)
tsne_out <- Rtsne(as.matrix(log2(t(lengthScaledTPM_qsmooth_highExp_remove_lowGenes_subTissue)+1)),perplexity = 40, check_duplicates = FALSE, theta=0.0 )
tsne_plot <- data.frame(tsne_out$Y)
tsne_plot$sample_accession <- colnames(lengthScaledTPM_qsmooth_highExp_remove_lowGenes_subTissue)
tsne_go <- tsne_plot %>% left_join(.,core_info)

Joining, by = "sample_accession"

set.seed(23235)
fit <- kmeans(tsne_go[,c('X1','X2')], 24)
cbind(tsne_go, fit$cluster) %>% group_by(fit$cluster,Tissue) %>% summarise(Count = n()) %>% mutate(freq = Count /sum(Count)) %>% summarise(Tissue_Count = n(), highest_freq=max(freq)) %>% mutate(mean(Tissue_Count), mean(highest_freq)) %>% select(`mean(Tissue_Count)`, `mean(highest_freq)`) %>% distinct()

OK, saving lengthScaledTPM_qsmooth_highExp_remove_lowGenes for downstream use:

save(lengthScaledTPM_qsmooth_highExp_remove_lowGenes, file='~/git/unified_gene_expression/data/lengthScaledTPM_processed.Rdata')